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Abstract. Adopting the theoretical framework for the generalized fishbonelike 
dispersion relation, an extended hybrid magnetohydrodynamics gyrokinetic simulation 
model has been derived analytically by taking into account both thermal ion 
compressibility and diamagnetic effects in addition to energetic particle kinetic 
behaviors. The extended model has been used for implementing an eXtended version 
of Hybrid Magnetohydrodynamics Gyrokinetic Code (XHMGC) to study thermal ion 
kinetic effects on Alfvenic modes driven by energetic particles, such as kinetic beta 
induced Alfven eigenmodes in tokamak fusion plasmas. 
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1. Introduction and motivation 

Nonlinear numerical simulations of magnetoliydrodynamics (MHD) and Alfven modes 
driven by energetic particles (EPs) mostly rely on hybrid MHD gyrokinetic codes, such 
as HMGC [I], M3D [2], and MEGA [3j. In the hybrid MHD gyrokinetic model, the 
thermal plasma component is described by MHD, while EP dynamics, in the so-called 
pressure coupling equation [2], is accounted for via the divergence of the EP pressure 
tensor, which is computed by solving the gyrokinetic equation with particle in cell (PIC) 
techniques. Kinetic treatments of the thermal plasma component are well known and 
generally implemented in (linear) spectral codes, such as NOVA-K O E] and MARS- 
K [7j. More recently, significant developments in gyrokinetic simulation codes, such as 
GTC [H] and GYRO [S], have also allowed investigating the kinetic effects of thermal 
plasma and EP dynamics on long wavelength electro-magnetic fluctuations, which 
were previously investigated only with codes based on the hybrid MHD-gyrokinetic 
approach PQ 121 El II]- In HMGC [T], the thermal plasma description is originally 
limited to the reduced MHD model [lO]. In the present work, our goal is to extend 
the hybrid MHD-gyrokinetic model implemented in HMGC to the low-frequency 
domain of the beta induced Alfven eigenmode (BAE) - shear Alfven wave (SAW) 
continuous spectrum [11], where the mode frequency can be generally comparable with 
thermal ion diamagnetic and/or transit frequencies, i.e. |ci;| ~ w^pj ~ uu. In this 
frequency range, where kinetic thermal ion (KTI) gap generally exists and influences 
plasma dynamics [TT], there is a continuous transition between various MHD and SAW 
fluctuation branches, as predicted theoretically [EllinilllllSlIISlIIZlIIHllIHlEO] and 
confirmed experimentally [2Il[221[231[21[2Sll2nil2Zll2Hl[2niEn]. Another notable feature 
of these low frequency fluctuations is that they may be resonantly excited by wave- 
particle interactions with EPs as well as thermal plasma particles, depending on the 
perpendicular wavelength jJOl El] . With the extended hybrid MHD gyrokinetic model 
discussed here, it will be possible to investigate various problems related with resonant 
excitation of Alfvenic and MHD fluctuations by EPs in the BAE-SAW continuous 
spectrum, consistent with gyrokinetic codes, e.g., GTC [32], in a common validity 
domain. Therefore, both the eXtended HMGC (XHMGC) and GTC codes can be 
verified using different models; yielding more detailed understanding of the underlying 
physics. In fact, theoretical and numerical work, presented in this article and partly 
developed within the framework of the SciDAC project on "Gyrokinetic Simulation 
of Energetic Particle Turbulence and Transport" (GSEP), was the prerequisite for 
successful verification of XHMGC predictions against analytic theories [33] as well as 
GTC numerical simulation results [SI EH] , reported recently. 

In this work, we extend the hybrid MHD-gyrokinetic model, derived originally 
in [2] for applications to numerical simulations of EP driven Alfven modes. The 
main differences with respect to the usual pressure coupling equation [2] are due 
to renormalization of the inertia term, to properly account for finite thermal ion 
diamagnetic effects, as well as to the gyrokinetic treatment of the thermal ion pressure 
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tensor, which allows us to properly handle wave-particle resonant interactions in the 
low frequency regime, where they can be of crucial importance for the analysis of linear 
and nonlinear behaviors of collisionless burning plasmas. The extended model has been 
developed assuming ideal Ohm's law as well as ignoring finite Larmor radius (FLR) 
effects in order to simplify the technical complications while still maintaining all essential 
physics ingredients [36]. In practice, maintaining the ideal MHD Ohm's law as limiting 
case implies assuming Tg <^ Tj and neglecting ion FLR effects, although finite magnetic 
drift orbit widths (FOW) are fully retained [37|. A more general approach without 
these simplifying assumptions will be developed in a separate work. For demonstrating 
the validity of the modified equations, we show that they are equivalent to the quasi- 
neutrality and vorticity equations derived in |36] for the frequency range from the 
kinetic ballooning mode (KBM) and BAE to the toroidal Alfven eigenmode (TAE). 
The XHMGC model equations in the linear limit are equivalent to the extended kinetic 
MHD used in spectral codes, such as NOVA-K O E] and MARS-K [Tj, but with EP 
dynamics treated non-perturbatively and on the same footing as the thermal plasma 
response (see Section [2] for more details). The possibility of investigating nonlinear 
dynamics, however, makes XHMGC more suitable to direct comparisons with M3D [2] 
or gyrokinetic codes [32l |38] in a common validity domain. 

The paper is organized as follows. In Section [2| the extended hybrid model 
equations are presented and discussed within the theoretical framework of [36]. In 
Section|3j we describe the numerical implementation of the extended model into HMGC, 
by adding both thermal ion compressibility and diamagnetic effects (of thermal ions as 
well as EPs) into MHD equations and a thermal ion population in the PIC module. 
In Section |4| possible applications and validity limits of XHMGC are discussed. A 
synthetic summary of current BAE numerical simulation results [3^ are also provided. 
Finally, conclusions and discussions are given in Section [5j 

2. Derivation of the extended hybrid model 

Reference [SB] presents a general theoretical framework for stability analyses of various 
modes and the respective governing equations. It shows that all modes of the shear 
Alfven branch having frequencies in the range between the thermal ion transit and 
Alfven frequency can be consistently described by one single general fishbone-like 
dispersion relation (GFLDR) [121 [HI EHl [El [H]. Reference [36] discusses various 
reduced equations governing the evolution of SAW fluctuations in burning plasmas, 
using the general approach of reference |39j. In this sense, reference [36] could seem not 
the optimal eliminate reference framework for further generalizing the HMGC hybrid 
model equations [11 HQ], which are to be used for nonlinear studies as well. However, 
the detailed analyses of reduced model equations, reported in |36j, on the basis of 
specializations of ordering of dimensionless parameters in the case of burning plasmas 
of fusion interest, starting from the somewhat different orderings of interest to space 
plasmas given in [52], allow to fully grasp the physics implications of the underlying 
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approximations. Moreover, on the basis of our discussions, it is straightforward to 
motivate the extension of the derived model equations to the nonhnear shown 
at the end of this section. 

Considering that the characteristic frequency, \u\, is much lower than the ion 
cyclotron frequency, \u!ci\, we may adopt the gyrokinetic theoretical approach and closely 
follow reference |39]. The low- frequency plasma oscillations can, thus, be described in 
terms of three fluctuating scalar fields: the scalar potential perturbation 6(j), the parallel 
(to b = B0/-B0, with Bo the equilibrium magnetic field) magnetic field perturbation SB\\ 
and the perturbed field Sip, which is related to the parallel vector potential fluctuation 
6A\\ by 

M|| = -i (^^) b ■ V6^. (1) 

The governing equations for describing the excitation of the shear Alfven frequency 
spectrum by energetic ions precession, precession-bounce and transit resonances in 
the range u^,pi ^ oJu < < uja, covering the entire frequency range from 
KBM/BAE [III [131 ini HI] to TAE |l2lll3llll], are generalized kinetic vorticity equation 
and quasi-neutrality condition, which can be written as followings, in the limit of 
vanishing FLR (see equation (16) and equation (17) in reference 
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(Aire \ An 

^ -^uu,sSK, ) + E TJr^^ X b ■ V(P,x + Ps\\)^M = 0, (2) 

-#) + Ee,(5ir,) = 0, (3) 

s=i 

where the non-adiabatic particle response, 5Ks., is obtained via the drift-kinetic equation 




[ojtrde - i{u} - uJd)]sSKs = M — ) QFqs 



{6(f) - 6ij) + 




(4) 



Here, angular brackets stand for velocity space integration, s denotes all particle 
species (e = bulk electrons, i = bulk ions, E = energetic particles), Cs and 
are the species electric charge and mass, Fq^ is the equilibrium distribution function 
(generally anisotropic), e = the energy per unit mass, QFqs = [uode + w*)sFos, 
'^*sFqs = ^cs^ X b)-VFos, k = — iV is the wave vector, Ucs = ^sBQ/msC is the cyclotron 
frequency, k± is the perpendicular wave vector, a;*ps = (k x b • V Pg) /nsmgOJcs is the 
diamagnetic frequency, Ps±_ and Pg^ are, respectively, the total perpendicular and parallel 
plasma pressures, Utr = v\\/qR is the transit frequency and Uds = ('^sc/es)(/i+t'j|/i?o)^K, 
with = k X b ■ K and k, = h ■ Vb. Note that the difference between Uds and 
^ds = {msc/es){fiflB + "Wy^KZ-Bo), with Qb = k x b ■ VBq/Bq, has been discussed 
in [361 ES] and, generally, must be handled properly; although, for many applications in 
low pressure {(3 = SttP/Bq <^ 1) plasmas, one can consider = Cods after solving for 
5B\\ from perpendicular pressure balance [39l|36|, as implicitly assumed in equations [21 13] 
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and |4} Note, also, that we have maintained the EP contribution to the divergence of 
the polarization current, which is represented by its leading term oc u^^pE in equation [2] 
This term is readily derived from the last term on the left hand side (LHS) of equation 



(13) in reference [3^ (see also Appendix A for further details) and was neglected in 



there due to the ordering I3e/ (^b ~ 'Te/tsd < 1, valid in a burning plasma dominated 
by fusion alpha particle self-heating. Here, (3e and /3b denote the beta values of EP and 
bulk plasma components (electrons and thermal ions), respectively, while te and tsd are 
the energy confinement time and EP slowing down time. More generally [T2l [T3l fHl 115] . 
the ordering Pe ~ Pb better represents nowadays magnetized plasmas of fusion interest 
and, thus, ue^jJ^pe ~ niU^pi, as assumed in equation |2| 

Equations [2] and [3| together with the drift-kinetic equation, equation |4} are the 
simplest yet relevant equations for analyzing the resonant excitations of SAW by EPs. 
Equation |2] demonstrates that both resonant as well as non-resonant responses due to 
the oc 6Ke term enter via the magnetic curvature drift coupling. In the high frequency 
case, oja > > ijJ*pi ^ oju, the thermal ion kinetic compression response 5Ki can be 
neglected. Thus, the quasi-neutrality condition, equation [3} reduces to the ideal MHD 
approximation, 50 ^ 5ip; i.e. 5E\\ ~ [SH]. Meanwhile, neglecting the oc oj^pi.oj^pE 
terms, equation |2] becomes equivalent to equation (3) in [2], i.e. the following pressure 
coupling equation in the hybrid MHD-gyrokinetic approach 

p,^ = -vn-(v.PE)x + i^; (5) 

where the subscript b denotes the bulk plasma (electrons and thermal ions), while pb 
and Vft are, respectively, bulk plasma mass density and fluid velocity. Here, the EP 
contribution to the perpendicular momentum change of the plasma has been neglected, 
due to He/ rib ^ \oj/uj*e\ [Ml 12], and thermal ion diamagnetic effects are consistently 
dropped since riEOJ^pE ~ ^iUj^pi. 

In order to extend the hybrid model to the low-frequency regime where uj ^ uu-, we 
need to include the effects of the thermal ion compressibility within the hybrid simulation 
scheme. That is, we need to include effects associated with the 5Ki terms in equations |2] 
and [3] First, in order to simplify the discussions, we formally assume Tg/Tj — )■ in 
the present work; the general case with finite Tg will be considered elsewhere. Thus, 
according to equation |3| we have (50 — (5^ ~ and the ideal MHD condition 5E\\ ~ 
remains valid. Next, we proceed to establish correspondences between the pressure 
coupling equation, equation [5| and the generahzed kinetic vorticity equation, equation [2] 

Applying the operator {d/dt)V ■ (Bo/-Bo)x to the linearized equation [s] and noting 
the quasi-neutrality condition V ■ J = 0, we readily derive 



(6) 



cdt " ' Bo ' cdt" ^ \BJ ' dt' \ Bl 




Extended hybrid model 



6 



Noting also the parallel Ampere's law along with V ■ 5 A = 0, 

47r(5J|| = -cV^M|j, (7) 

and equation [l| term (i) can be seen to correspond to the field line bending term; 
i.e. the first term in equation |2} Term (ii), on the contrary, does not have any 
direct correspondence in equation [2j This term is the usual kink drive and it was 
dropped in the analysis of [36j, focusing on drift Alfven fiuctuations with high mode 
numbers, for it is formally of 0{l/n), with n the toroidal mode number. However, 
as noted in equation (Al) of [3^, term (ii) is readily recovered in a form that can be 
straightforwardly reduced to that reported here. Meanwhile, from the linearized Ohm's 
law 

5Ex + -5vb X Bo = 0, (8) 

c 



and SE± = — V_l50, term (iii) corresponds to the second term in equation |2j with 
the oc (jj^pi.uj^pE terms neglected. To establish correspondences between the pressure 
responses in equations |2] and [6| we first denote = Pg + Pj. It can then be shown 
(see appendix Appendix B ) that term (iv) corresponds to the thermal ion and electron 
contributions to the last term on the LHS of equation [2| when kinetic compression effects 
of the background thermal plasma are neglected. 

Finally, let us discuss term (v), due to EP pressure perturbation, which can be 
expressed as (see Appendix C). 

of ■ ( Bo ) = iS""'*^-' + 

Meanwhile, noting the definition of 5Ks [391 EE], the 5Ks term in equation [2] can be 
shown to be related with the pressure perturbations as (see Appendix D ) 

\k^c^ I kjcBQ 



4:71 

+ b) ■ (VPo.x + VPos\\)nM. (10) 

Note that the 2nd term in the right hand side (RHS) disappears in the ideal MHD 
5(j) ^ 6tp limit. Equation 10, thus, clearly demonstrate that the SKg contribution in 
equation [2| combined with the 3rd term on the RHS of equation 10 (or the last term on 
the LHS in equation |2|, has the same form of equation [o] and recovers the total pressure 
response of term (iv) in equation [6] for SP±i = SP^ = SPi. In other words, the SKg 
term in equation |2] corresponds to the kinetic compressibility component of the pressure 
perturbations. 

Summarizing the above discussions, it is clear that, in order to include effects 
due to finite thermal ion compressibility and diamagnetic drift as well as the finite 
EP contribution to the divergence of the polarization current, the pressure coupling 
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equation in the MHD-gyrokinetic approach, equation |5} has to be modified such that its 
perpendicular components are given by equation A.9 of [Appendix A , which we rewrite 
here for the reader's convenience 



Pb 



dt 



+ V;, • V + 



b X VPqe^ 



- V^Pe - (V ■ P,)^ - (V ■ P, 



+ 



J X B 



c 



Here, v;, = b x VPoi±/{pbijJci) + ^^b, ^^b = {c/Bq)6E x b and the "unshifted" pressure 
tensors Pe and Pi need to be calculated from solutions of the gyrokinetic equations 
as specified in [Appendix A while Pg is consistently neglected in the present approach, 
assuming Tg/Tj 0. Reminding the concluding remark of Appendix A, this equation 
readily reduces to the well-known pressure coupling equation |2} in the limit where 
thermal ion diamagnetic effects and EP contribution to the divergence of the polarization 
current are neglected. 

As anticipated above, in the present work, we followed reference ^36j , since that has 
a detailed discussion of vahdity hmits of different reduced models of the whole vorticity 
and quasi-neutrality equations, derived for fusion applications and following the trace 
of reference 
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includes equilibrium parallel current effects, as discussed 
earlier in this section and in [36] (Appendix). This simple remark readily follows from 
the discussion presented in [31] as well as the modified momentum balance equation 



implemented in XHMGC, i.e. equation 11 itself. The present model is valid in the 
nonlinear case too, as shown by the simple derivation provided in Appendix A and by 
the following discussion. This is deduced easily from direct inspection of equation (5) 
in reference [1^. That equation clearly shows that, for the small FLR limit considered 
in HMGC [H 130], the nonlinear terms, treated explicitly, are those that are coming 
from convective E x B nonlinearity and from the Maxwell stress nonlinearity, when 
the thermal ion response is taken in the fiuid limit, both of which are readily obtained 
from equation 11 upon application of the operator dtV ■ (Bq/Bq)x, as it was done 
for equation [5] earlier in the section. Other nonlinear dynamics, which are implicitly 
included in (V ■ Pi) and (V ■ Pe)± terms, are fully retained via equation 11 Thus, 



the back reaction of zonal structures onto SAW fiuctuations is fully accounted for, 
i.e. that of zonal fiows (ZFs) and fields as well as radial modulation of equilibrium 
profiles [01 [m, HE] which also enter via the diamagnetic terms in equation 11 computed 
on the whole (slowly evolving) thermal ion and EP pressure profile, obtained from 
the respective toroidally and poloidally averaged distribution functions. This choice is 
consistent with known approaches to nonlinear MHD equations, accounting for finite 
diamagnetic drift corrections [ITJ HHl SOI EOl EI] • 

Thus, the approximations involved with the extended implementation within 
XHMGC on the basis of equation 11 consist of neglecting FLR, assuming electron as 
a massless fiuid, considering Te/Tj — )■ (such that parallel Ohm's law is recovered 
in the ideal MHD limit) and accounting for Reynolds stress in the thermal ion fiuid 
limit. The possible further extension of the present model to include finite Tg/Tj and 
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generalizing the parallel Ohm's law, while maintaining other simplifying assumptions, 
is straightforward on the basis of the present discussion and will be reported in a 
separate work [52]. Here we note that the present extended hybrid model, based on 



equation 11, with clearly formulated assumptions that limit its applicability, includes 
very rich physics; e.g. it is capable to correctly evaluate the renormalized inertia for ZFs, 
for which the trapped thermal ion dynamics is of crucial importance, and to account for 
geodesic acoustic mode (GAM) kinetic response, including Landau damping. 

So far, XHMGC has been used for moderate EP drive [331 [Ml EHl [37] , where the EP 
diamagnetic correction to the divergence of the polarization current can be neglected, 
as argued in [3S]- Actually, in the studies reported in [33], thermal ion diamagnetic 
contribution to the polarization current is also neglected, since the case of uniform 
thermal ion pressure profiles is investigated in there for facilitating comparisons of 



numerical simulation results with analytic theory predictions (see also section |4]). 
3. Numerical implementation 

HMGC [1] is used for investigating linear and nonlinear properties of moderate toroidal 
number (n) shear Alfven modes in tokamaks. It solves the coupled set of O(e^) reduced- 
MHD equations [10] for the electromagnetic fields and the gyro-center Vlasov equation 
for a population of energetic ions, where large aspect ratio is assumed, i.e. e = a/Ro <^ 1, 
with a and Rq the tokamak minor and major radius, respectively . Energetic particles 
contribute to the dynamic evolution of the wave fields via the pressure tensor term in the 
MHD equations, as described by the pressure coupling equation [2] . This code allows us 
to describe both self-consistent mode structures in toroidal equilibria and EP dynamics, 
as well as to get a deeper insight into how the Alfvenic modes affect the confinement of 
such particles. 

The extended model, described in section [2| has been implemented into the 
extended version of HMGC (XHMGC). Following the general procedure, described in 



references [HHQ], for the formal manipulation of equation 11, the relevant equations for 
the MHD solver are in terms of the poloidal magnetic field stream function \Ef and U, 
which is proportional to the scalar potential $ and defined as U = —c^/Bq, can be 
written in the following form in the cylindrical coordinate system [R, Z, ip): 

- ^ -V* X . VC/ + + + 0(.S.S,), (12, 

J D 2 dU 
^ [Dt ^ RqOZ 



R^ dPn, 




+ 



i^RoUcio dZ ' R^iUcEo dZ 
_ ^ ( R' dP,,^ ^ R' dP,E^ \ 

\Rlojcio dZ RloJcEo dZ ) ^ 
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B 



VPoE± , _ VPoj± 

Vv9 X h X 



VViU 



Rq 




(13) 

where we have maintained the same notation of reference [1] and exphcitly show the 
additional terms that have been added to implement the extended XHMGC model. 
Thus, 

. i?2 D d R^ 
dt^R 



P 



rCr. 



Dt 



-VU X V(/? ■ V, 



the Grad-Shafranov operator A* is defined by 
~ ORRdR^ dZ^' 

Bq is the vacuum magnetic field on the magnetic axis at = -Rq [H ^ao = ^iBo/ {mic), 
<^cEo = ^eBqI {niEc) and the subscript _L denotes components perpendicular to ip. In 
the above equations, v_|_ is the E x B fluid velocity, p is the bulk plasma mass density, 
Pe is the scalar pressure of bulk electrons, II^; and Ilj are, respectively, the pressure- 



(14) 



tensor of the EP and thermal ions, computed with the definitions of equations |A.3 
given in 

CO 



and |A.7 



Appendix A and, r] is the resistivity and c is the speed of light. 

hmited to the MHD 



These O(e^) equations have been first derived in reference 
description of the thermal plasma, while the inclusion of energetic particle dynamics 
has been discussed in references. P, HQ]. Here, in the proposed further extension of the 
numerical model, thermal ion dynamics as well as diamagnetic effects are also taken 
into account, according to equation 11, derived in Appendix A and section |2j At the 
leading order in e, O(e^), the reduced-MHD equations describe the thermal plasma in 
the cylindrical approximation. Toroidal geometry enters the equations as corrections at 
the next order in the inverse aspect ratio. 

In order to close equations 12 and 13, the EP and thermal ion pressure tensor 
components can be obtained by directly calculating the appropriate velocity moments 
of the distribution function for the particle population interacting with the perturbed 
electromagnetic field. As discussed in section [2| we initially assume the Tg/Tj — )■ limit 
for the sake of simplicity, i.e. Pe — )■ 0. Meanwhile, with cold electron assumption and 
ignoring thermal ion finite Larmor radius (FLR), ideal MHD parallel Ohm's law can be 
readily recovered. 

As to numerical formulation, the equations of motion in gyro-center coordinates 
for thermal ions are in the same form, mutatis mutandis, as those reported in [T] for 

X Please, note the difference between tlie present notation, wliere Bq stands for the on axis equihbrium 
magnetic field, and that used in section II, where Bq generally denoted the (spatially dependent) 
equilibrium magnetic field. 
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EPs. In the gyrocenter-coordinate system Z = (R, M, V, 6), where R is the gyrocenter 
position, M is the conserved magnetic moment, V is the parallel speed and 6 is the 
gyrophase, the equations of motion take the form 



dR 

~dF 



dM 
dt 
dV 
~dt 



Vh + 



h X V0 



V 



-b X Van 



+ 



M V 



V + 



m. 



b X V In 5, 



+ 



-b-{ 



V 



Van X V(/)} - 



V0 + 



M . 



-Van 



X Vln5 



b- VlnB. 



(15) 



Here, the subscript s denotes either EP or thermal ion species and, using the same 
notations as in [T], fig = esBQ/mgC is the corresponding cyclotron frequency. The 
fluctuating potential ay is related to the poloidal magnetic field stream function \E' 
through the relationship ay = [eg/ c){Rq/ R)'^ . The parallel electric field term in the 
equation for V has been suppressed, neglecting, thus, small resistive corrections to the 
ideal-MHD parallel Ohm's law. Meanwhile, the pressure tensor can be written, in terms 
of the gyrocenter coordinates, as 

4 / dZDz^^2Fs(t,%M,V) 



n.(t,x) 



mi 



m,. 



I + bb 



y2 _ ^ 



(5(x - R) 



(16) 



where I is the unit tensor, lij = 6ij, Fs(t, R, M, "K) is the gyrocenter distribution 
function and D^^^z is the Jacobian of the transformation from canonical to gyrocenter 
coordinates. The distribution function Fs satisfies the Vlasov equation 




V + 



dV d ' 
~dtdV 



0, 



(17) 



where dH/dt and dV /dt are given by equation 15 In the numerical implimentation 



of XHMGC, equations [15] and [17] can be readily solved as a full-F simulation. On the 
other hand, a 5f algorithm [101 ISSl EH ES] is also implemented in order to minimize the 
discrete particle noise. The latter is recommended as far as 5f <C F, the former when 
5f ^ F. 



4. Applications 

In general, the extended version of HMGC can have two species of kinetic particles. On 
one hand, one can use XHMGC for investigating thermal ion kinetic effects on Alfvenic 
modes driven by EP. On the other hand, it may be interesting to use XHMGC as a tool 
to simulate two coexisting EP species, generated e.g. by both ion cyclotron resonance 
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heating (ICRH) and neutral beam injection (NBI) heating, in order to study hnear 
excitation of Alfvenic fluctuations and Energetic Particle Modes (EPM) [15] , as well as 
the interplay between the respective nonlinear physics controlled by the different heating 
sources [56] . 

HMGC has been extensively used in pj and |57j to investigate the linear physics 
(damping and EP drive mechnisms), and in [30] and [58] to analyze the nonlinear 
dynamics of EPM. XHMGC has been verified against those previous findings and can 
recover numerical simulation results in the above studies. Furthermore, by accounting 
for the kinetic thermal ion effects, XHMGC shows the existence of Kinetic BAE (KBAE) 
which can be seen as radially trapped eigenstates due to discretization of BAE-SAW 
continuum by FLR/FOW effects, as well as KBAE resonantly excited by wave-particle 
interactions with EPs [33} [37]. 

As an example to demonstrate the capability of XHMGC, we briefly report 
simulation results of KBAE, which is discussed in detail in [33]. The results show that 
a fully kinetic treatment of thermal ions is necessary for a proper description of the low 
frequency Alfvenic fluctuation spectrum. By including thermal ion compressibility, our 
numerical simulations do show the existence of a finite-frequency BAE accumulation 
point in the SAW continuum, which was demonstrated analytically and numerically 
using MHD codes [59]. Meanwhile, when effects due to finite ion drift orbit width 
(FOW) are included, our simulations clearly demonstrate that the BAE-SAW continuum 
becomes discretized; yielding a series of discrete kinetic eigenmodes with small frequency 
separation [33l [60] . In figure [l| we have plotted the BAE accumulation frequencies in 
the fluid limit which are defined as ubae = (l^ti{^ + Te/TiY/"^ with Tg/Tj ^ in the 
current case, as well as eigenmode frequencies obtained from simulation results. The 
analytically predicted KBAE frequencies [33] are in good agreement with observations 
from numerical simulations. The results also indicate that FOW kinetic effects increase 
with the toroidal mode number, as expected [331 [60] . 

On the other hand, our simulations also show that KBAE can be driven by EPs. 
In figure |2| we can see that the frequencies scale properly with the KBAE frequencies; 
and the growth rates decrease with the thermal ion temperature due to the stronger ion 
Landau damping and/or the weaker EP drive due to the increased frequency mismatch 
between mode and characteristic EP frequencies. In the absence of thermal ion kinetic 
effects, the excited modes may be identified as energetic particle mode (EPM); which 
requires sufficiently strong drive to overcome the SAW continuum damping. Including 
the thermal ion kinetic effects not only introduce a finite kinetic thermal ion frequency 
gap at the BAE accumulation frequency but also discretize the BAE-SAW continuum. 
In that case, the continuum damping is greatly reduced or nullified, and the discrete 
KBAE's are more readily excited by the EP drive. 
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Figure 1. Real frequency comparison between simulations by "antenna" excitations 
and the theoretical accumulation point frequencies. Simulations refer to an equilibrium 
magnetic field charatcrized by shifted circular magnetic siirfaccs with inverse aspect 
ratio a/Ro = 0.1 and the q-profile, in the cylindrical approximation, given by 
q{l) = q{0) + [q{l) - q{0)]r^ , where r is normalized to a, q{0) = 2.7 and q{l) = 3.9. A 
is eigenfrequency from simulations for n = 1, + is eigenfrequency from simulations for 
n = 3, the black dashed line is cobae the accumulation point frequency. 




Figure 2. The real frequency u) and growth rate 7 for the n=3 mode versus 
different thermal ion pressure parameters for Pi = 0.0072,0.0128,0.02, and with a 
fixed value of (3e = 0.009. "*" is the mode real frequency of simulation results by EP 
excitations; "A" is the KBAE frequencies by antenna excitations; solid line denotes 
the theoretical BAE accumulation point frequency; "+" is the growth rate by EP 
excitation simulations. 
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5. Conclusions and discussions 

In the present work, we have employed the theoretical framework (generalized 
kinetic vorticity and quasi-neutrality equations) of the generalized linear fishbone 
dispersion relation and derived an extended hybrid MHD-gyrokinetic simulation model 
applicable to the low-frequency regime, where effects of thermal ion compressibility and 
diamagnetic drifts play significant roles in the dynamics of Alfven waves and energetic 
particles in tokamak plasmas. The extended simulation model has been implemented 
into an eXtended version of HMGC (XHMGC). Initial simulations of XHMGC have 
discovered the existence of KBAE discretized by the thermal ion FOW effects, which 
are absent in conventional MHD codes. Simulations also demonstrate that KBAE can 
be readily excited by EPs. In the current model, we have taken Te/Tj — )■ and neglected 
finite Larmor radius effects in order to simplify the presentation and focus on the most 
important qualitative new physics connected with implementation of the thermal ion 
compressibility. In addition, XHMGC is limited to circular shifted magnetic surfaces 
equilibria, with relatively large aspect ratio; XHMGC includes kinetic effects related to 
both bulk and fast ions; however, it is typically used for retaining only the perturbed 



pressure for two EP species; XHMGC doesn't include rotation (see Appendix A), while it 
retains the perturbed electrostatic potential. These additional effects will be considered 
in future works. 

More recently, the electromagnetic formulation ^ of global gyrokinetic particle 
simulation in toroidal geometry has been implemented in GTC [32] • In such a code, 
ions are treated by the gyrokinetic equation, while electrons are simulated using an 
improved fluid-kinetic electron model [H]. In [31], the connection between the extended 
hybrid MHD-gyrokinetic model and gyrokinetic simulation model has been verified in 
the drift kinetic limit as well as ignoring the terms on the order of 0((e/g)^). Instead 
of directly calculating the pressure tensor, lower moments of the kinetic equation have 
been calculated, i.e. the perturbed density and parallel current. Using charge neutrality 
condition, it can be demonstrated that the combination of the perturbed density and 
parallel current contribution is totally equivalent to the pressure tensor in equation [TT] 
Therefore, both GTC and XHMGC can be verified using different models in a common 
validity regime; yielding more detailed understanding of the underlying physics. 

On the other hand, kinetic compressibility is also included in (linear) spectral codes, 
such as NOVA-K [SIE] and MARS-K [7]. While in NOVA-K only perturbative treatment 
of the kinetic effects is considered and continuum damping is not taken into account, as 
in MARS-K, since both are eigenvalue codes, the spectral approach allows the study of 
the linear stability of all eigenmodes in a general equilibria; meanwhile, kinetic effects 
are generally related to bulk plasmas only, although the inclusion of fast ions is quite 
straightforward. As to other hybrid MHD gyrokinetic codes, M3D [2] is based on the 
pressure coupling equation; MEGA [3l H] uses a hybrid model for MHD and energetic 
particles, where the effect of the energetic ions on the MHD fluid is taken into account 
in the MHD momentum equation through the energetic ion current. The diamagnetic 



Extended hybrid model 



14 



drift effect is evaluated in the MHD equations by adding the diamagnetic advection 
term to the equation of motion [HI SHI SH |50l [61]. At present, XHMGC can handle 
two species kinetic particles self-consistently, but is limited to circular shifted magnetic 
flux surfaces equilibria with vanishing bulk plasma equilibrium pressure. Meanwhile, a 
new version of the code with general equilibria is being developed, with the capability 
of solving fully compressible gyrokinetic particle response. 
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Appendix A. Simple derivation of model equations 

Adopting a multi-fluid moment description of plasma dynamics, the force balance 
equation can be written as 

V6 + Pi? I ^ + Vij ■ V I Ve 

= -VPe-V-P,-V-Ps + ^^^ . (A.l) 

c 

Here, pb and pe are bulk plasma and EP mass densities, v;, = b x VPoi±/{pb<^ci) + S^b, 
= hx'VPoE±/{PE^cE)+S^b+^UE\\ and 6\b = (c/i?o)5Ex b from equation (8), having 
omitted terms that are 0{co^:pe/i^ce) or higher with respect to the RHS. Furthermore, 
thermal ion and EP pressure tensors on the RHS have to be interpreted as usual, i.e. 
with the conventional fluid velocity shift in the definition 



Psij =^s J d\{Vi - Usi){Vj - Usj)fs , (A. 2) 

with fs the particle distribution function and Usi = J dwvifs/ris. When the pressure 
tensor is computed form the particle distribution function within the gyrokinetic 
description, some subtleties are connected with the ordering Usi/vts ~ Pls/L in the 
plane orthogonal to b, with p^g the Larmor radius of the s-species, Vts its thermal speed 
and L the characteristic equilibrium radial scale-length. Thus, in the drift-kinetic limit 
used in this work, = Ps_lI + {Ps\\ — Ps_L)bb, with I the unit diagonal tensor and 

2 

Ps± = m,j d^j'^fs , (A.3) 

Ps\\ ='^s j dw{v\\ - Us\\offs ■ (A.4) 
Note the difference between Us\\q = J d-vv\\Fos/nos, used here, and Us\\ = J dwufs/ns, 



used in equation A. 2 Fqs being the (slowly evolving) equilibrium particle distribution 



function. In equation |A.1[ we assumed that only EPs can carry significant parallel fluid 
velocity. 
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The perpendicular component of equation A.l can be further simphfied, by noting 
that Vb = hK,{l + 0(e^)), with e = a/i?o and a and Rq the tokamak minor and major 
radii, and that we are using the optimal ordering |w| ~ |w*pi| <C |w*p£;|- This allows us 
to rewrite 



PE 



dt 



+ ^E-^\^E 



b X VPm± 



i^cE 



(A.5) 



where we have dropped the oc {dt + ■ V)^^E± terms, for they are 0{u/co^:pE) 
and, similarly, the oc M_B||b ■ Vv£;_l term, since it is ©[(Tj/Tg;)^/^] - or, equivalently, 
0([?2£;/nfe)^/^] - for w ~ UJ^:pi ~ cuu, at shorter wavelength or higher frequency, this term 
would be negligible anjrway with respect to the thermal ion inertia response, as negligible 



would be diamagnetic responses of both EPs and thermal ions. So, equation A.5 well 
describes the physics we want to incorporate in the present analysis. Recalling the 
definition of P^, we also have 

- (V ■ P,)^ = -V±Ps± - k{Ps\i - Ps±) (A.6) 

Thus, the second term on the RHS of equation A.5 can be combined with the oc Pe\\ 



term on the RHS of equation |A.6[ computed for EPs, and actually be reabsorbed into 
that (up to the relevant order), provided that the pressure tensor is reinterpreted as 
P5x)bb, with the "unshifted" expression 



P. = Ps±l + (Psl 

p.l 



(A.7) 



replacing the usual definition given in equation A. 4 With this convention on the 



pressure tensor, the perpendicular components of equation A.l can be rewritten as 



Pb 



' d_ 
dt 



+ Vb ■ V V5 



b X VPoE± 

i^cE 



Vv 



E± 



- V±Pe - (V ■ P.)^ - (V 




(A. 



Actually, equation A.8 can be reduced further when residual terms that are 0{u}^pe/^ce) 
or higher with respect to the RHS are omitted, as noted below equation |A.1[ In fact, 
one readily obtains 



pb 



'd_ 
dt 



+ Vfe ■ V + 



b X VPoEL 

^cE 



V 



5\h 



VxPe - (V ■ P.)^ - (V ■ Ve)^ + 



J X B 



(A.9) 



This equation readily reduces to the well-known pressure coupling equation [2], in the 
limit where thermal ion diamagnetic effects and EP contribution to the divergence of 
the polarization current are neglected. 

Appendix B. Study of term (iv) in equation [6] 



In the low-/3 approximation (Vlni?o — k) 
' b \ 2b X K 



V X 



(B.l) 
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Meanwhile, in the incompressible limit, 

d 



dt 



5Pb + (5vb ■ VPc 



06 



0, 



where 



Svh± = c 



Bo X V±5^ 



R2 



Then, with equations B.l B.2 and B.3 

b 



a V X 



5n 



2b X K 
2b X K 



9t 



■ V( 



Bp X Vj 

-^0 



■ VPr 



06 J 



-"0 



So 

X b ■ VPobflJ 



where fi^ = k x b ■ k. 



(B.2) 
(B.3) 



(B.4) 



Appendix C. Study of term (v) in equation [6] 

Assuming 

(5Pe = hhSPEW + (I - hh)6PE±, 

we can show 

V-6Pe = {6Pei\ - SPE±)ihV -h + K) 

+ hVii{6PEii-6PE±) + \/6PE±, 
where k = h ■ Vb. Thus, 

b X V ■ = b X V6Pe± + {6Pe\\ - SPE±)h x k. 

Now 



V ■ (— X V6Pe±) 

-Do 



2b X K 

Sn 



£±5 



and 



V ■ [(5P^,|| - 5Pe±)^ xk]^ ■ V(5Ps|| - (5Pijx) 

-Do -Do 



(C.l) 

(C.2) 
(C.3) 
(C.4) 

(C.5) 



In equation C.5[ we have used the large aspect ratio assumption, e ^ 1, consistent with 



the reduced MHD description used in HMGC [HllQ]. Combining equations C.3 to C.5 
we obtain 



^ / bx (V-(5Pe) 
dt '[ Bo 



Bn 



n.{6PEii + sPE±] 



(C.6) 
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Here, we assume the definition of [32] 
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(D.i; 



de CO 

where 6fs is the fiuctuating particle distribution function, , on the RHS, we have dropped 
all terms oc dF^s/dfi, for they generate contributions of higher order in what follows [36j. 
Thus, in our treatment, Fqs is generally anisotropic, although terms oc dFos/dfi do not 
appear explicitly. One then finds 



Aires 
Hc2 



Os 



de 



5<p 



II 



+ 



J^{k^Ps))S^, (D.2) 



III 



with < ■ ■ ■ > denoting velocity integration and Jq is the zero order Bessel function. 
For term (I), we obtain 



where 



and 



6P. 



Attu 

knC 



m 



vl5fs 



2 



(D.3) 
(D.4) 
(D.5) 



are, respectively, perturbed perpendicular and parallel pressure. 
Meanwhile, for \k±p\ <^ 1, term (III) can be written as 

„,2 



(///) 



47re 



Anel 



klc ' B 



^^^^ dF^ 

klmoC^ de 



esB 
{6(f) - dtp) 



k X b) ■ VF( 



Os 



5tp 



An 



(D.6) 



Combining equations D.3, D.6 and term (II) with the first term on the RHS of 



equation D.6 we obtain 
/ Aire 



\ ky 



A-nel 
^ kgmsC^ 
Att 



dF 



UUds 



Os 



de 



+ tj7^^k{\^ X b) ■ (VPos± + VPo.||)5V^. (D.7) 
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